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The nonlinear saturation of the r-mode instability and its effects on the spin evolution of Low 
Mass X-ray Binaries (LMXBs) are modeled using the triplet of modes at the lowest parametric 
instability threshold. We solve numerically the coupled equations for the three mode amplitudes in 
conjunction with the spin and temperature evolution equations. We observe that very quickly the 
mode amplitudes settle into quasi-stationary states that change slowly as the temperature and spin 
of the star evolve. Once these states are reached, the mode amplitudes can be found algebraically and 
the system of equations is reduced from eight to two equations: spin and temperature evolution. The 
evolution of the neutron star angular velocity and temperature follow easily calculated trajectories 
along these sequences of quasi-stationary states. The outcome depends on whether or not the 
star will reach thermal equilibrium, where the viscous heating by the three modes is equal to 
the neutrino cooling (H = C curve). If, when the r-mode becomes unstable, the star spins at a 
frequency below the maximum of the H = C curve, then it will reach a state of thermal equilibrium. 
It can then either (1) undergo a cyclic evolution with a small cycle size with a frequency change 
of at most 10%, (2) evolve toward a full equilibrium state in which the accretion torque balances 
the gravitational radiation emission, or (3) enter a thermogravitational runaway on a very long 
timescale of » 10 6 years. If the star does not reach a state of thermal equilibrium, then a faster 
thermal runaway (timescale of ~ 100 years) occurs and the r-mode amplitude increases above the 
second parametric instability threshold. Following this evolution requires more inertial modes to 
be included. The sources of damping considered are shear viscosity, hyperon bulk viscosity and 
viscosity within the core-crust boundary layer. We vary proprieties of the star such as the hyperon 
superfiuid transition temperature T c , the fraction of the star that is above the threshold for direct 
URCA reactions, and slippage factor, and map the different scenarios we obtain to ranges of these 
parameters. We focus on T c > 5 x 10 9 K where nonlinear effects are important. Wagoner [1[ has 
shown that a very low r-mode amplitude arises at smaller T c . For all our bounded evolutions the 
r-mode amplitude remains small ~ 10" 5 . The spin frequency of accreting neutron stars is limited by 
boundary layer viscosity to ^ max ~ 800-ffz[S lla /(Mi.47?6)] 4 ^ 11 r g ~ 2/ ' 11 • Fast rotators are allowed for 
[S ns /(Mi. 4 iie)] 4/U T g - 2/11 ~ 1 and we find that in this case the r-mode instability would be active 
for about 1 in 1000 LMXBs and that only the gravitational waves from LMXBs in the local group 
of galaxies could be detected by advanced LIGO interferometers. 

PACS numbers: 04.40.Dg, 04.30.Db, 97.10.Sj, 97.60.Jd 



I. INTRODUCTION 

R-modes are oscillations in rotating fluids that are due 
to the Coriolis effect. They are subject to the classical 
Chandrashekar-Friedman-Shutz (CFS) instability @, Q, 
which is driven by the gravitational radiation backreac- 
tion force. Andersson [j] and Friedman and Morsink Q 
showed that, in the absence of fluid dissipation, r-modes 
are linearly unstable at all rotation rates. However, in 
real stars there is a competition between internal viscous 
dissipation and gravitational driving Q that depends on 
the angular velocity Q and temperature T of the star. 
Above a critical curve in the fl—T plane the n = 3, m = 2 
mode, referred to as 'the r-mode' in this work, becomes 
unstable. At first, an unstable r-mode grows exponen- 
tially, but soon it may enter a regime where other in- 
ertial modes that couple to the r-mode become excited 
and nonlinear effects become important. Roughly speak- 
ing, nonlinear effects first become significant as the am- 
plitude passes its first parametric instability threshold, 
which is very low (~ 10~ 5 ). Modeling and understand- 
ing the nonlinear effects is crucial in determining (1) the 



final saturation amplitude of the r-mode and (2) the lim- 
iting spin frequency that neutron stars can achieve. The 
r-mode amplitude and the duration of the instability are 
among the main factors that determine whether the as- 
sociated gravitational radiation could be detectable by 
laser interferometers on Earth. 

The r-mode instability has been proposed as an expla- 
nation for the sub-breakup spin rates of both Low Mass 
X-ray Binaries (LMXBs) @, i] and young, hot neutron 
stars d, Q . The idea that gravitational radiation could 
balance accretion was proposed independently by Bild- 
sten and Andersson et al. Q. Cook, Shapiro and 
Teukolsky [HEH model the recycling of pulsars to mil- 
lisecond periods via accretion from a Keplerian disk onto 
a bare neutron star with M = 1.4M Q when O = 0. De- 
pending on the equation of state they found that spin 
frequencies of between ss 670 Hz and 1600 Hz could be 
achieved before mass shedding or radial instability set 
in (these calculations predated the realization that the 
r-mode instability could limit the spin frequency). For 
comparison, the highest observed spin rate of millisec- 
ond pulsars is 716 Hz for PSR J1748-2446ad [jj, 0- 
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PSR B1937+21, which was discovered in 1982, was the 
previous fastest known radio pulsar with a spin rate of 
642 Hz [3]; that this "speed" record stood for 24 years 
suggests that neutron stars rotating this fast are rare. 
Moreover, based on a Bayesian statistical analysis of the 
spin frequencies of the 11 nuclear-powered millisecond 
pulsars whose spin periods are known from burst oscil- 
lations, Chakrabarty et al. [151 ] claimed a cutoff limit of 
^roax = 760 Hz (95% confidence); A more recent analy- 
sis, which added two more pulsars to the sample, found 
iW = 730 Hz 0. 

At first sight, one might conclude that mass shedding 
or radial instability sets ^ max , and that it is just above 
the record v = 716 Hz determined for PSR J1748-2446ad. 
However, the nuclear equations of state consistent with 
this picture all have rather large radii ~ 16—17 km 
for non-rotating 1.4 M Q models; see Table 1 in Cook et 
al. [l0(. For these equations of state, the r-mode insta- 
bility should lead to zAnax somewhat below 716 Hz; see 
Eq. (|33|) in Sec. [V] below. Thus, the r-mode instability 
may prevent recycling by accretion from reaching mass 
shedding or radial instability. In other words, the de- 
tection of the 716 Hz rotator is consistent with accretion 
spin-up mitigated by the r-mode instability only for equa- 
tions of state for which mass shedding or radial instability 
would permit even faster rotation. Ultimately, this may 
be turned into useful constraints on nuclear equations of 
state. However, at present the uncertainty in the physics 
of internal dissipation is a significant hindrance in estab- 
lishing such constraints. 

Since a physical model to follow the nonlinear phase of 
the evolution was initially unavailable, Owen et al. [ItJ 
proposed a simple one-mode evolution model in which 
they assumed that nonlinear hydrodynamics effects satu- 
rate the r-mode amplitude at some arbitrarily fixed value. 
According to their model, once this maximum allowed 
amplitude is achieved, the r-mode amplitude remains 
constant and the star spins down at this fixed ampli- 
tude (see Eqs. (3.16) and (3.17) in Ref. [l7j]). They used 
this model to study the impact of the r-mode instability 
on the spin evolution of young hot neutron stars assum- 
ing normal matter. In their calculation they include the 
effects of shear viscosity and n-p-e bulk viscosity. They 
found that the star would cool to approximately 10 9 K 
and spin down from a frequency close to the Kepler fre- 
quency to about 100 Hz in a period of ~ 1 yr [17] . 

Most subsequent investigations that did not perform 
direct hydrodynamic simulations used the one-amplitude 
model of Ref. [l7| for studying the r-mode instability. 
Levin [18[ used this model to study the limiting effects of 
the r-mode instability on the spin evolution of LMXBs, 
assuming an r-mode saturation amplitude of ~ 1; he 
adopted a modified shear viscosity to match the maxi- 
mum LMXB spin frequency of 330 Hz known in 1999. 
Levin found that the neutron star followed a cyclic evo- 
lution in the f2 — T phase plane. The star spins up for 
several million years until it crosses the r-mode stability 
curve, whereupon the r-mode becomes unstable and the 



star is viscously heated for a fraction of a year until the 
r-mode reaches its saturation amplitude (~ 1). At this 
point the spin and r-mode amplitude evolution equations 
are changed, following the prescription of Ref. [l7| to en- 
sure constant amplitude. The star then spins down by 
emitting gravitational radiation for another fraction of a 
year until it crosses the r-mode stability curve again and 
the instability shuts off. The time period during which 
the r-mode is unstable was found to be about 10~ 6 times 
shorter than the spin-up time, and Levin concluded that 
it is unlikely that any neutron stars in LMXBs in our 
galaxy are currently spinning down and emitting gravita- 
tional radiation. However, following work by Arras et al. 
[l9| showing that nonlinear effects become significant at 
small r-mode amplitude, Heyl (20| varied the saturation 
amplitude, and found that the duration of the spin-down 
depends sensitively on it. He predicted that the unstable 
phase could be as much as 30% of the cyclic evolution for 
an r-mode saturation amplitude of a w 10~ 5 , and that 
this would make some of the fastest spinning LMXBs in 
our galaxy detectable by interferometers on Earth. 

Jones [2l| and Lindblom and Owen [22j pointed out 
that if the star contains exotic particles such as hyperons 
(massive nucleons where an up or down quark is replaced 
with a strange quark), internal processes could lead to a 
very high coefficient of bulk viscosity in the cores of neu- 
tron stars. While this additional high viscosity coefficient 
could eliminate the instability altogether in newly born 
neutron stars [2l|, [2^, [H, [24[ , Nayyar and Owen [24[ pro- 
posed that it would enhance the probability of detection 
of gravitational radiation from LMXBs by blocking the 
thermal runaway. 

The cyclic evolution found by Levin [l8| and gener- 
alized by Heyl [20j | arises when shear or boundary layer 
viscosity dominates the r-mode dissipation. In the evo- 
lutionary picture of Nayyar and Owen [24|, the r-mode 
first becomes unstable at a temperature where shear 
and boundary layer viscosity dominate, but the result- 
ing thermal runaway halts once hyperon bulk viscosity 
becomes dominant. The key feature behind the runaway 
is that shear and boundary layer viscosities both decrease 
with increasing temperature, so the instability speeds up 
as the star grows hotter. However, if the bulk viscosity 
is sufficiently large the star can cross the r-mode stabil- 
ity curve at a point where the viscosity is an increas- 
ing function of temperature. Such scenarios were stud- 
ied by Wagoner 1] for hyperon bulk viscosty with low 
hyperon superfluid transition temperature; similar evo- 
lution was found for strange stars by Andersson, Jones 
and Kokkotas [25j. In this picture, the star evolves near 
the r-mode stability curve until an equilibrium between 
accretion spin-up and gravitational radiation spin-down 
is achieved. The value of the r-mode amplitude remains 
below the lowest instability threshold found by Brink et 
al. [261 l27l |28| for modes with n < 30, and hence in this 
regime nonlinear effects may not play a role. 

Schenk et al. [2^| developed a formalism to study the 
nonlinear interaction of the r-mode with other inertial 
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modes. They assumed a small r-mode amplitude and 
treated the oscillations of the modes with weakly nonlin- 
ear perturbation theory via three-mode couplings. This 
assumption was tested by Arras et al. 19] and Brink 
et al. [H, [27], [28| . Arras et al. proposed that a turbu- 
lent cascade would develop in the strong driving regime. 
They estimated that r-mode amplitude was small and 
could have values between 10 _1 — 10~ 4 . Brink et al. 
modeled the star as incompressible and calculated the 
coupling coefficients analytically. They computed the in- 
teraction of about 5000 modes via approximatively 1.3 
million couplings of the 10 9 possible couplings among 
the modes with n < 30. The couplings were restricted to 
mode triplets with a fractional detuning 5ui/(2£l) < 0.002 
since near-resonances promote modal excitation at very 
small amplitudes. Brink et al. showed that the nonlinear 
evolution saturates at a very small amplitude, generally 
comparable to the lowest parametric instability thresh- 
old that controls the initiation of energy sharing among 
the sea of inertial modes. However, Brink et al. did not 
model accretion spin-up or neutrino cooling in their cal- 
culation and only included minimal dissipation via shear 
viscosity. 

In this paper we begin a more complete study of the 
saturation of the r-mode instability including accretion 
spin up and neutrino cooling. We use a simple model in 
which we parameterize uncertain properties of the star 
such as the rate at which it cools via neutrino emission 
and the rate at which the energy in inertial modes dis- 
sipates via boundary layer effects [30( and bulk viscos- 
ity. In order to exhibit the variety of possible nonlinear 
behaviors, we explore a range of models with different 
neutrino cooling and viscous heating coefficients by vary- 
ing the free parameters of our model. In particular, we 
vary: (1) the slippage factor 5 ns , which regulates the 
boundary layer viscosity, between and 1 (see for exam- 
ple [U [ill, HH for some models of the interaction between 
the oscillating fluid core and an elastic crust) ; (2) the 
fraction of the star that is above the density threshold for 
direct URCA reactions /dUj which is taken to be between 
(0% of the star cools via direct URCA) and 1 (100% of 
the star is subjected to direct URCA reactions), and in 
general depends on the equation of state used; and (3) the 
hyperon superfluidity temperature T c , which is believed 
to be between 10 9 — 10 10 K (We use a single, effective T c 
rather than modelling its spatial variation.) We focus on 
T c > 5 x 10 9 K for which nonlinear effects are important. 
For low T c < 3 x 10 9 K, Wagoner [l[ showed that the 
evolution reaches a steady state at amplitudes below the 
lowest parametric instability threshold found by Brink 
et al. [28]. It is important to note that all our evolu- 
tions start on the part of the r-mode stability curve that 
decreases with temperature and that the bulk viscosity 
does not play a role in any of our bound evolutions. 

We include three modes: the r-mode at n — 3 and the 
two inertial modes at n = 13 and n = 14 that become 
unstable at the lowest parametric instability threshold 
found by Brink et al. [281 ] . We evolve the coupled equa- 



tions for the three-mode system numerically in conjunc- 
tion with the spin and temperature evolution equations. 
The lowest parametric instability threshold provides a 
physical cutoff for the r-mode amplitude. In all cases we 
investigate, the growth of the r-mode is initially halted 
by energy transfer to the two daughter modes. We ob- 
serve that the mode amplitudes settle into a series of 
quasi-stationary states within a period of a few years af- 
ter the spin frequency of the star has increased above the 
r-mode stability curve. These quasi-stationary states are 
algebraic solutions of the three-mode amplitude equa- 
tions (see Eqs. ©) and change slowly as the spin and 
the temperature of the star evolve. Using these solutions 
for the mode amplitudes, one can reduce the eight evo- 
lution equations (six for the real and imaginary parts of 
the mode amplitudes, which are complex [291 ]: one for the 
spin, and one for the temperature) to two equations gov- 
erning the rotational frequency and the temperature of 
the star. Our work can be regarded as a minimal physical 
model for modeling amplitude saturation realistically. 

The outcome of the evolution is crucially dependent 
on whether the star can reach a state of thermal equilib- 
rium. This can be predicted by finding the curve where 
the viscous heating by the three modes balances the neu- 
trino cooling, referred to below as the Heating = Cooling 
(H = C) curve. The H = C curve can be calculated prior 
to carrying out an evolution using the quasi-stationary 
solutions for the mode amplitudes. If the spin frequency 
of the star upon becoming unstable is below the peak 
of the H = C curve, then the star will reach a state of 
thermal equilibrium. When such a state is reached we 
find several possible scenarios. The star can: (1) un- 
dergo a cyclic evolution; (2) reach a true equilibrium in 
which the accretion torque is balanced by the rate of loss 
of angular momentum via gravitational radiation; or (3) 
evolve in thermal equilibrium until it reaches the peak of 
the H — C curve, which occurs on a timescale of about 
10 6 yr, and subsequently enter a regime of thermal run- 
away. On the other hand, if the star cannot find a state 
of thermal equilibrium, then it enters a regime of ther- 
mogravitational runaway within a few hundred years of 
crossing the r-mode stability curve. When this happens, 
the r-mode amplitude increases beyond the second para- 
metric instability, and more inertial modes would need 
to be included to correctly model the nonlinear effects. 
This will be done in a later paper. 

This paper focuses on showing how nonlinear mode 
couplings affect the evolution of the temperature and spin 
frequency of a neutron star once it becomes prone to the 
r-mode CFS instability. We do this in the context of three 
mode coupling, which may be sufficient for large enough 
dissipation. To illustrate the types of behavior that arise, 
we adopt a very specific model in which the mode fre- 
quencies and couplings are computed for an incompress- 
ible star, modes damp via shear viscosity, boundary layer 
viscosity and hyperon bulk viscosity, and the star cools 
via a mixture of fast and slow processes. This model in- 
volves several parameters that are uncertain, and we vary 
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these to find 'phase diagrams' in which different generic 
types of behavior are expected. Moreover, the model it- 
self is simplified: (1) A more realistic treatment of the 
modes could include buoyant forces, and also mixtures of 
superfluids or of superfluid and normal fluid in different 
regions. (2) Dissipation rates, particularly from bulk vis- 
cosity, depend on the composition of high density nuclear 
matter, which could differ from what we assume. 



Nevertheless, although the quantitative details may 
differ from what we compute, we believe that many fea- 
tures of our calculations ought to be robust. More sophis- 
ticated treatment of the modes of the star will still find a 
dense set of modes confined to a relatively small range of 
frequencies. Most importantly, this set will exhibit nu- 
merous three mode resonances, which is the prerequisite 
for strong nonlinear effects at small mode amplitudes. 
Thus, whenever the unstable r-mode can pass its lowest 
parametric instability threshold, it must start exciting its 
daughters. Whether or not that occurs depends on the 
temperature dependence of the dissipation rate of the r- 
mode; for the models considered here, where bulk viscos- 
ity is relatively unimportant, soon after the star becomes 
unstable its r-mode amplitude passes its first paramet- 
ric instability threshold. Once that happens, the generic 
types of behavior we find - cycles, steady states, slow 
and fast runaway - ought to follow suit. The details of 
when different behaviors arise will depend on the precise 
features of the stellar model, but the principles we out- 
line here (parametric instability, quasisteady evolution, 
competition between heating and cooling) ought to ap- 
ply quite generally. 

In Sec. |TT] we describe the evolution equations of the 
three modes, the angular frequency and the temperature 
of the neutron star. We first show how the equations of 
motion for the modes of Schenk et al. couple to the rota- 
tional frequency of the star in the limit of slow rotation. 
We then give a short review of the parametric instability 
threshold and the quasi-stationary solutions of the three- 
mode system. The thermal and spin evolution of the star 
is discussed next. This is followed by a description of 
the driving and damping rates used. Sec. |TTT] provides 
an overview of the results, which includes a discussion 
of each evolution scenario and of the initial conditions 
and input physics that lead to each scenario. Sec. IIV Al 
discusses cyclic evolution in more detail. An evolution 
that leads to an equilibrium steady state is presented 
next in Sec. IIVBI The two types of thermal runaway are 
then discussed in Sec. IIV CI The prospects for detecting 
gravitational radiation for the evolutions in which the 
three-mode system correctly models the nonlinear effects 
are considered in Sec. |V] We summarize the results in 
the conclusion. Appendix A sketches a derivation of the 
equations of motion for the three modes and Appendix 
B contains a stability analysis of the evolution equations 
around the thermal equilibrium state. 



II. EVOLUTION EQUATIONS 

A. Three mode system: coupling to uniform 
rotation 

In this section we review the equations of motion for 
the three-mode system in the limit of slow rotation. In 
terms of rotational phase r for the time variable with 
dr = Q dt Eq. (2.49) of Schenk et al. [29[ can be rewritten 
as 
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Here the scaled frequency Qj is defined to be ujj = w^/ft, 
the dissipation rates of the daughter modes are and 7 7 , 
"f a is the sum of the driving and damping rates of the r- 
mode 7 Q = jgr — lav, and the dimensionless coupling is 
k = k/ (MR 2 il 2 ). These amplitude variables are complex 
and can be written in terms of the variables of Ref. [29[ 
as Cj(t) = y/fl(i)cj(t) (see Appendix A for a derivation 
of Eqs. [I]). The index j loops over the three modes j = 
a, /3, 7, where a labels the r-mode or parent mode and 
and 7 label the two daughter modes in the mode triplet. 

When the daughter mode amplitudes are much smaller 
than that of the parent mode, one can approximate the 
parent mode amplitude as constant. Under this assump- 
tion one performs a linear stability analysis on Eqs. ([1]) 
and finds the r-mode amplitude when the two daughter 
modes become unstable (see Eqs. (B5-B7) of Ref. [28[ 
for a full derivation). This amplitude is the parametric 
instability threshold 
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where the fractional detuning is Su = Co a — Cop — <1> 7 . 
Thorough explorations of the phase space of damped 
three-mode systems were performed by Dimant 3_4| and 
Wersinger et al. [35l ]. 

For the three modes at the lowest parametric instabil- 
ity threshold, u a « 0.66, wp w 0.44, Co 1 w 0.22, h w 0.19 
and \Suj\ « 3.82 x 10 -6 . Note that Co is twice the w of 
Brink et al. [H, H3, HI • Here /3 labels the mode with 
n = 13, m = —3 and 7 labels the n = 14, m = 1 mode. 
The amplitude the r-mode has to reach before exciting 
these two daughter modes is \C a \ w 1.5 x lO^-s/ft [28| . 

We next rescale the rotational phase t by the fractional 
detuning as f — t\5£j\ and the mode amplitudes by 



|C a |o 
|C 7 |o 



\Sid\yTT c \5Q\Vn. 



\8Co\\/Vr c 



4Ky / 0) ct CJ 7 ' 



(3) 



5 



which for the r-mode is, up to a factor of 
the no-damping limit of the parametric instability thresh- 
old below which no oscillations will occur. The coupled 
equations become 



dC a 
~dY 

dC 
df 

dC-y 

~df~ 



Following Eq (K39-K42) of Schenk et al. [2g| the physical 
angular momentum of the perturbation can be written as 

ttJphys = C B°A J d 3 X P [(Cl X C B ) ■ (A >< U) (») 
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with Cj = Cj-/|(7j|o and jj = 7j/n c being the newly 
rescaled amplitudes and dissipation/driving rates, re- 
spectively. 



Quasi- Stationary Solution 



In terms of amplitudes and phase variables Cj 
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Eqs. (0| can be rewritten as 
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Since the eigenvectors £4 oc e mA ^ the cross-terms will 
vanish for modes with different magnetic quantum num- 
bers to as J e i( - mA - mB ^d(j) = for toa 7^ m B . Eq. © 
can be re-written for our triplet of modes as 

J phys = MR 2 {k aa \C a \ 2 + k p0 \C p \ 2 + fc 77 |C 7 | 2 ), (9) 

where k aa is defined as 

K a = — |jp | rf 3 ^p[(r2 x £*) ■ (O x £ a ) - i«aCi • (^ x £«)] 

( 10 ) 

and similarly for fc/3/3 and A: 77 . In terms of the scaled vari- 
ables Cj = Cj/\Cj\ (with \Cj\ defined in Eq. ©) the 
angular momentum of the perturbation can be written 
as 
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where we have defined the relative phase difference as 
<t> = 4> a — 4>f3 — 4>j. These equations have the stationary 
solution 
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We chose the same normalization for the eigenfuctions 
as Refs. 0, [H, H3, [H, so that at unit amplitude all 
, modes have the same energy e a — MR 2 ft 2 . The energy 
of a mode a is E a = MR 2 n 2 \c a \ 2 = MR 2 VL\C a \ 2 . The 
rotating frame energ y is the same as the canonical energy 
and physical energy [29] . The canonical angular momen- 
tum and the canonical energy of the perturbation satisfy 
the general relation E c = —(uj /m)J c 

Angular momentum is gained because of accretion and 
lost via gravitational waves emission 

d.J 



1 + 

7/3 ' 



tan^ 

- 77 - 7a 



— = 2 lGR J c rmodc + M VGMR, (12) 
dt 

where J crmo dc = -(m a /uj a )e a \c a \ 2 = -3Afi? 2 Sl|c Q | 2 = 
-3MR 2 \C a \ 2 . Eq. ((12]) can be rewritten in terms of the 
scaled variables C« as 



n\5Q\ 



Note that in the limit in which 7^+7 7 >> j a the station- 
ary solution for the r-mode amplitude \C a \ is the same 
as the parametric instability threshold. 



B. Temperature and Spin Evolution 

The spin evolution equation is obtained from conser- 
vation of total angular momentum J, where 



dJ _ 67GR MR 2 n c \5Cj\ 



df Q {Ak) 2 uif)C 



C a \ 2 + 
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(13) 



J 1*1 H~ t/phys- 
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Thermal energy conservation gives the temperature evo- 
lution equation 

dT 

C(T)— = Y, 2E 'jl3+ K nMc 2 -L V {T), (14) 

3 

= 2MR 2 n( lav \C a \ 2 + l0 \C p \ 2 
+ ^\C 7 \ 2 ) + K n Mc 2 - L V (T). 
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The three terms on the right hand side of the equa- 
tion represent viscous heating, nuclear heating and neu- 
trino cooling. The specific heat is taken to be C(T) w 
1.5 x 10 38 T 8 ergKT 1 , where T = T 8 x 10 8 K. Nu- 
clear heating occurs because of pycnonuclear reactions 
and neutron emission in the inner crust 36]. At large 
accretion rates such as that of the brightest LMXBs of 
M w lO -8 M0/yr, the accreted helium and hydrogen 
burns stably and most of the heat released in the crust is 
conducted into the core of the neutron star, where neu- 
trino emission is assumed to regulate the temperature of 
the star [36|, [37| ■ The nuclear heating constant is taken 
to be K n « 1 x 10" 3 Following Ref. [J, we take the 
neutrino luminosity to be 

L v = LdijT® Rd\j(T/Tp) + L m \jT g R m \j(T/T p ) (15) 



where the constants for the modified and direct URCA re- 
actions are defined by L m \i = l.Ox 10 32 erg sec -1 , L^u = 
/du x 10 8 L m u [H, |3{|, and the electron-ion, neutron- 
neutron neutrino bremsstrahlung and Cooper pairing of 



neutrons are given by £ e -i = 9.1 x 10 29 erg sec 1 
L n _ n « 0.01L mU , L Cp = 8.9 x 10 31 erg sec" 1 [40]. The 
fraction of the star /du that is above the density thresh- 
old for direct URCA reactions is in general dependent on 
the equation of state [4l[ and in this work we treat /du 
a free parameter with values between and 1. 

The proton superfluid reduction factors for the modi- 
fied and direct URCA reactions are taken from Ref. [39[ 
(see Eqs. (32) and (51) in Ref. [H): 

5.5 



R dV (T/T p ) = 0.2312 + V(0-76880) 2 + (0.1438v) 2 (16) 

x exp ^3.427 - v/(3.427) 2 + v 2 ^j , 

(T/T p ) = (0.2414 + ^/(O^Se) 2 + (0.1318u) 2 ) 7 , 

x exp (5.339 - V(5-339) 2 + (2v) 



R 



where the dimensionless gap amplitude v for the singlet 
type superfluidity is given by 





1 1.456 - 0.157^/^ + 1.764^ J . (17) 



Similar to Ref. [lj], we use T p — 5.0 x 10 9 K. In terms of 
the scaled variables Eq. (|14[) becomes 



C(T)— = p \ uj alav \C a \ 2 +u 0l p\Cp\ 2 18 

ar [fin) u} a u)[}U)~ 



+&J% |C 7 I') + 



K n Mc 2 -L v (T) 



C. Temperature and Spin Evolution with the 
Mode Amplitudes in Quasi-Stationary States 

Assuming that the amplitudes evolve through a series 
of spin- and temperature-dependent steady states, i.e., 



dCi/df ~ 0, the spin and thermal evolution equations 
can be rewritten by taking J rs ID, and using Eqs. © in 
Eq. HT 



dft 
~df 
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kaa I 1 



1 



(19) 



M VGMR 
r%MR 2 i£l\6u 



where / = I /(AIR 2 ). The thermal evolution of the sys- 
tem is given by 



df (4re) 2 a; a ,w 1 aw 7 Q,\Su\ \ la 

1 \ K n M<?-L v (T) 



up (20) 



1 



tan J 



n c n\scj\ 



By setting the right hand side of the above equation to 
zero, one can find the Heating = Cooling (H = C) curve. 
Below, we find that Eqs. (TTS"]) - (|2T))) describe the evolu- 
tion very well throughout the unstable regime. These 
equations are a minimal physical model for the effects of 
nonlinear coupling on r-mode evolution. 



D. Sources of Driving and Dissipation 

The damping mechanisms are shear viscosity, bound- 
ary layer viscosity and hyperon bulk viscosity; for modes 
j = a, (3, 7 we write 

7, „(n, T) = 7j sh (T) + 7j bl (fi, T) + jj hb (fi, T). (21) 

The r-mode is driven by gravitational radiation and 
damped by these dissipation mechanisms, while the pair 
of daughter modes (n = 13, m = —3 labeled as (3 and 
n = 14, to = 1 labeled as 7) is affected only by the vis- 
cous damping. Brink et al. [26|, [2?], [28| determined that 
this pair of modes is excited at the lowest parametric 
instability threshold. Their model uses the Bryan [42j 
modes of an incompressible star, which has the advan- 
tage that the mode eigenfrequencies (and eigenfunctions) 
are known analytically. This enables them to find near 
resonances efficiently We are using their results, but we 
include more realistic effects such as bulk viscosity, whose 
effect vanishes in the incompressible limit (ri — > 00 in 
Eq. j29l)) 

For our benchmark calculations, we adopt the neutron 
star model of Owen et al. Ref. [6j (n = 1 polytrope, 
M = 1.4M©, n c = 8.4 x 10 3 rad sec" 1 and R = 12.53 
km) and use their gravitational driving rate and shear 
viscous damping rate for the r-modc 



sec 



O 6 
3.26 

7ash(T) ~ —2, 

T s h J-s 



(22) 



7 



where r s h = 2.56 x 10 6 sec. (In Sec. [V] we consider ap- 
proximate scalings with M and R.) 

The damping rate due to shear viscosity for the two 
daughter modes is calculated using the Bryan modes for 
a star with the same mass and radius 

7 /3 sh (T) ~ 3.48 x lO^sec" 1 ^, (23) 



7 7 sh(T) 



L 8 

4.52 x 1(T 4 scc" 1 -^-. 



The geometric contribution "fsh/v of the individual modes 
increases significantly with the degree n of the mode scal- 
ing approximatively like n 3 for large n (see Eq. (29) of 
Brink et al. [27j for an analytic fit to the shear damping 
rates computed for the 5,000 modes in their network), 
and hence the inertial modes with n = 13 and n = 14 
have shear damping rates about three orders of magni- 
tude larger than that of the r-mode. 

The damping due to boundary layer viscosity is calcu- 
lated using Eq. (4) of Ref. [H, 



7abl(7 1 , fl) 

7 /3 bi(7 1 , n) 
7 7 bi(T, n) 



0.009 sec" 1 S 2 



0.028 sec -1 Si 



(24) 



0.021 sec" 1 ^ v 
is 



Analogous to Wagoner [l|, we allow the slippage fac- 
tor Sps to vary. The slippage factor is defined by Refs. 
[E WM to be S* s = [2Sl + 5 s 2 )/3, with S n being the 
fractional difference in velocit y o f the normal fluid be- 
tween the crust and the core [3l|] and S s the fractional 
degree of pinning of the vortices in the crust [45] . Note 
that 7^ bi and 7 7 y are both greater than 2 x j a w and 
can easily be comparable to 7gr in the unstable regime. 

The damping rate due to bulk viscosity produced by 
out-of-equilibrium hyperon reactions for the r-mode is 
found by fitting the results of Nayyar and Owen [24| . 
This rate is taken to have a form similar to that taken 
by Wagoner [l[ 



la hb = fl 



hb - 



' 1 + {uMt(T))2 ' 
and for the daughter modes 



7/3 hb = f\ 



hb " 



t pT{T)C 



•0 



'l + (%Or(T)) 2 ' 
and similarly foe 7 7 hb- The relaxation timescale 



(25) 



(26) 



(27) 



Rhb(T/T c ) 

The reduction factor is taken to be the product of two 
single-particle reduction factors [H, [24| 

a 5/4 + b l/2 
Rhb single (T /T c ) = 



w here a = 1 + 0.3118?/ 2 , 6 = 1 + 2.556y 2 and y = 
^1.0 - T/T C (1.456 - 0.157v/T c /T + 1.764T C /T). The 
constants t\ w 10~ 4 sec and tp n ~ 0.00058 sec are found 
by fitting the results of Ref. 24] . The factor /hb allows 
for physical uncertainties; we take /hb = 1 throughout 
the body of the paper since T c , which enters jj hb ex- 
ponentially, is also uncertain. For the daughter modes, 
the dissipation energy due to bulk viscosity is calculated 
using the modes for the incompressible star. In the slow 
rotation limit, it is given to leading order in T^ 2 by 



exp (0.5068 - v/0.5068 2 + y : 
(28) 




p 



(29) 



This approximation was proposed by Cutler and Lind- 
blom [43| and adopted by Kokkotas and Stergioulas [44| 
for the r-mode and by Brink et al. [13] for the inertial 
modes. The adiabatic index Ti is regarded as a parame- 
ter; we use fi w 2. The damping rate is 



E 



(30) 



where e = MR 2 il 2 is the mode's energy in the rotat- 
ing frame at unit amplitude and j = (3,^. Using this 
procedure, we calculate 



'0/3 ' 



1.4 x 10 5 sec 
: 1.0 x 10~ 5 



(31) 



sec. 



III. SUMMARY OF RESULTS 

Fig. DJa) shows possible evolutionary trajectories of a 
neutron star in the angular velocity-temperature — Tg 
plane, where T = T§ x 10 8 K is the core temperature, and 
il = fl/Cl c = fl j \JnGp with p the mean density of the 
neutron star. Fig. [ljb) displays the regions in /du — S ns 
in which the trajectories occur. Here /du represents the 
fraction of the star that is above the density threshold 
for direct URCA reactions and S ns is the slippage factor 
that reduces the relative motion between the crust and 
the core taking into account the elasticity of the crust 
[3TI ] . The stability regions are shown at fixed hyperon 
superfluidity temperature, T c = 5.0 x 10 9 K. The initial 
part of the evolution is similar in all scenarios and can 
be divided into phases. 

Phase 0. Spin up below the r-mode stability curve at 
Tg, = Tgi n such that nuclear heating balances neutrino 
cooling. 

Phase 1. Linear regime. The r-mode amplitude grows 
exponentially. The phase ends when the r-mode reaches 
the parametric instability. 

Phase 2. The triplet coupling leads to quasi-steady 
mode amplitudes. The star is secularly heated at 
approximately constant Q because of viscous dissipation 
in all three modes. 

Phase 3. Several trajectories are possible depending on 
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FIG. 1: (a)Typical trajectories for the four observed evolu- 
tion scenarios are shown in the Cl - Ts phase space, where 
SI = 0,/Q.c. The dashed lines (H = C curves) represent 
the points in the Q — Ts phase space where the dissipative 
effects of the heating from the three-modes exactly compen- 
sate the neutrino cooling for the given set of parameters (S ns , 
fdu, T c , ...) of each evolution. (b)The corresponding sta- 
bility regions for which these scenarios occur are plotted at 
fixed hyperon superfluidity temperature T c — 5.0 x 10 9 K, 
while varying /du and S ns - The position of the initial angu- 
lar velocity and temperature (Qi n , Is i n ) with respect to the 
maximum of this curve determines the stability of the evo- 
lution. (I) f2i n > Ojj^cmax- Trajectory Ri. Fast Runaway 
Region. After the r-mode becomes unstable the star heats 
up, does not find a thermal equilibrium state and continues 
heating up until a thermogravitational runaway occurs. (II) 
Qin < £Ih=c max- The evolutions are either stable or, if there 
is a runaway, it occurs on timescales comparable to the ac- 
cretion timescale. The possible trajectories are (l)Trajectory 
C. Cycle Region. (2) Trajectories Si and S2. Steady State 
Region. (3) Trajectory i?2- Slow Runaway Region. 



how the previous phase ends. 

a. Fast Runaway. The star fails to reach thermal 
equilibrium when the trajectory passes over the peak of 
the Heating = Cooling (H = C) curve. This leads to 
rapid runaway. The daughter modes damp eventually 
as bulk viscosity becomes important, and the r-mode 
grows exponentially until the trajectory hits the r-mode 
stability curve again. This scenario ends as predicted 
by Nayyar and Owen (24|. However, the r-mode passes 
its second parametric instability threshold soon after it 
starts growing again. This requires the inclusion of more 
modes to follow the evolution, which is the subject of 
future work. 

b. The star reaches thermal equilibrium. There are then 
three possibilities: 

(i) Cycle. The star cools and spins down slowly, 
descending the H = C curve until it crosses the r-mode 
stability curve again. At this point the instability shuts 
off. The star cools back to T§i n at constant (l and 
then the cycle repeats itself. At T c = 5.0 x 10 9 K this 
scenario occurs for values of S ns < 0.50 and large enough 
values of /du- However, if T c is larger, the cycle region 
in the jdirSns phase space increases dramatically (see 
Fig- Eta))- Note that our cycles are different from those 
obtained by Levin [l8| in that the spin-down phase 
does not start when the r-mode amplitude saturates 
(or in our case when it reaches the parametric insta- 
bility threshold), but rather when the system reaches 
thermal equilibrium. The r-mode amplitude does not 
grow significantly above its first parametric instability 
threshold, remaining close to ~ 10 5 and so the part 
of the cycle in which the r-mode is unstable also lasts 
longer than in Ref. [18j. Also, our cycles are narrow. 
During spin-down the temperature changes by less than 
20 % and changes by less than 10% of the initial value. 
(See Sec. [2] for a detailed example.) 

(ii) Steady State. For small S ns and large enough 
/du (/du > 5 x 10- 5 , S ns < 0.04; see Fig. Hp))) the 
star evolves towards an f2 equilibrium. The trajectory 
either ascends or descends the H = C curve (spins 
up and heats or spins down and cools). The evolution 
stops when the accretion torque equals the gravitational 
radiation emission. 

(iii) Slow Runaway. For small S ns and very small /du 
(S ns < 0.03, /du < 5 x 10~ 5 ) the star ascends the H = C 
curve until the peak is overcome and subsequently a 
runaway occurs. The daughter modes eventually damp 
and the r-mode grows exponentially until it crosses its 
second parametric instability threshold and more modes 
need to be included. 

Bulk viscosity only affects the runaway evolutions; the 
cyclic and steady state evolutions found here would be 
the same if there were no hyperon bulk viscosity. For 
large T c ~ 10 10 , or for models with no hyperons at all, 
there would be no runaway region (See Fig. [DJa) for an 
/du — S ns scenario space with a larger T c — 6.5 x 10 9 K 
where the fast runaway region has shrunk dramatically 
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FIG. 2: Two cyclic trajectories in the Q, — Ts plane are dis- 
played for a star with T c = 5.0 x 10 9 K and (a) /au = 0.15 
and Sn 8 = 0.10, and (b) / du = 0.142 and S ns = 0.35, which is 
close to the border between the stable and unstable region (see 
Fig.fjjb)). The thick solid line labeled as the Heating = Cool- 
ing (H = C) curve is the locus of points in this phase space 
where the neutrino cooling is equal to the viscous heating due 
to the unstable modes. The other solid line representing the 
r-mode stability curve is defined by setting the gravitational 
driving rate equal to the viscous damping rate. The part of 
the curve that decreases with Ts is dominated by boundary 
layer and shear viscosity, while the part of the curve that has 
a positive slope is dominated by hyperon bulk viscosity. In 
portion a\ — > fei of the trajectory the star heats up at con- 
stant H. Part 6i — > c\ represents the spin down stage, which 
occurs when the viscous heating is equal to the neutrino cool- 
ing. c\ — > di shows the star cooling back to the initial T$. 
Segment di — > a± displays the accretional spin-up of the star 
back to the r-mode stability curve. The cycle a-2, — > di pro- 
ceeds in the same way. This cycle is close to the peak of the 
H = C curve. Configurations above this peak will run away. 



and the slow runaway region has disappeared.) 



IV. POSSIBLE EVOLUTION SCENARIOS 

In this section we examine examples of the different 
types of evolution in more detail. We assume M = 
Kr 8 Af Q /yr and T c = 5.0 x 10 9 K. 



A. Cyclic Evolution 

In this sub-section we present the features of typical 
cyclic trajectories of neutron stars in the angular velocity 
temperature plane in more detail. We focus on two cases: 
(CI) S ns = 0.10 and / dU = 0.15 and (C2) 5 ns = 0.35 
and /du = 0.142. In this scenario the 3-mode system is 
sufficient to model the nonlinear effects and successfully 
stops the thermal runaway. The numerical evolution is 
started once the star reaches the r-mode stability curve. 
The initial temperature of the star is at the point where 
nuclear heating equals neutrino cooling in Eq. (|18[) that 
is approximately Tg,i n ~ 3.29 for both cases. The initial 
51 is the angular velocity that corresponds to this tem- 
perature on the r-mode stability curve, which differs for 
the different S ns (O in = 0.183 for d and Sl in = 0.288 for 
C 2 ). 

Figs. (Ha) and (b) display the cyclic evolution for tra- 
jectories Ci and C*2 of Fig. [TJb). In leg ai — > b\ of the 
trajectory the r-mode and, once the r-mode amplitude 
increases above the first parametric instability thresh- 
old, the two daughter modes it excites, viscously heat 
up the star until point b\ when the neutrino cooling bal- 
ances the viscous dissipation. This part of the evolu- 
tion occurs at constant angular velocity over a period 
of ^hcat-up ~ 100 yr and a total temperature change 
(AT) ai _ bl w 0.80 (w 24% of T 8in ). The points where 
the viscous heating compensates the neutrino cooling are 
represented by the Heating = Cooling (H = C) curve. 
This is determined by setting Eq. [TH] to zero and using 
the quasi-stationary solutions given by Eq. ^ for the 
three modes on the right hand side. The star continues 
to evolve on the H = C curve for part b\ — > c\ of the 
trajectory as it spins down and cools down back to the r- 
mode stability curve. This spin-down stage lasts a time 
ispin-downbi-ci ~ 23,000 yr that is much longer than 
the heat-up period. This timescale is very sensitive to 
changes in the slippage factor and can reach 10 6 yr for 
smaller values of ,5ns that are close to boundary of the 
steady state region. The cycle is very narrow in angular 
velocity with a total angular velocity change of less than 
4%, (Af2)(, 1 _ Cl w 0.0066. The temperature also changes 
by only about 2%, (AT 8 )/ Jl _ Cl « 0.08 in this spin-down 
period. Segment ci — > d\ represents the cooling of the 
star to the initial temperature on a timescale of ~ 2, 000 
yr. In part d± — > a\ the star spins up by accretion at 
constant temperature back to the original crossing point 
on the r-mode stability curve. This last part of the tra- 
jectory is the longest-lasting one, taking w 200, 000 yr 
at our chosen M of 10 -8 M Q yr -1 . The cycle in Fig. 
[2Jb) proceeds in a similar fashion. It is important to note 



that this configuration is close to the border between the 
"FAST RUNAWAY" and "CYCLE" regions and there- 
fore close to the peak of the H — C curve. Configura- 
tions above this peak (e.g., with the same /du and higher 
Sns) w ih g° through a fast runaway. 

Fig. OHa) shows the evolution of the three modes in 
the first few years after the star first reaches the r- 
mode stability curve. In this region the r-mode is un- 
stable and initially grows exponentially. Once it has in- 
creased above the first parametric instability threshold 
the daughter modes are excited. The oscillations of the 
three modes display some of the typical dynamics of a 
driven three-mode system. When the r-mode transfers 
energy to the daughter modes they increase exponen- 
tially while the r-mode decreases. Similarly, when daugh- 
ter modes decrease the r-mode increases. The viscosity 
damps the oscillations and the r-mode amplitude settles 
at a value close to the parametric instability threshold. 
Fig. E^b) displays the evolution of the r-mode ampli- 
tude divided by the parametric instability threshold on 
a longer timescale. It can be seen that the r-mode never 
grows significantly beyond this first threshold. Fig. [3^c) 
shows the evolution of the parametric instability thresh- 
old as a function of time. The threshold increases as the 
temperature increases and the star is viscously heated by 
the three modes. When the star spins down in thermal 
equilibrium, the threshold decreases to a value close to 
its initial value. 




1 .05 
1 

0.95 
0.9 
_ 0.85 
a 0.8 

_o 

-a 0.75 
_o 

0.7 

0.65 
0.6 
0.55 

n 




0.6 



0.7 



0.N 



0.9 



i [x 10 /8w] 



B. Steady State Evolution 



This sub-section focuses on evolutions that lead to a 
steady equilibrium state in which the rate of accretion 
of angular momentum is balanced by the rate of loss 
via gravitational radiation emission. This scenario is re- 
stricted to stars with small slippage factor (5 ns < 0.04, 
see Fig. [ljb)) and boundary layer viscosity. A typical 
trajectory of a star that reaches such an equilibrium is 
shown in Fig. 2] As always, we start the evolution at the 
point on the r-mode stability curve at which the nuclear 
heating balances neutrino cooling. Above the r-mode sta- 
bility curve the gravitational driving rate is greater than 
the viscous damping rate and the r-mode grows exponen- 
tially until nonlinear effects become important. In this 
case, as in the cyclic evolution, the triplet of modes at 
the lowest parametric instability threshold is sufficient to 
stop the thermal runaway. The r-mode remains close to 
the first instability threshold for the length of the evo- 
lution and after a few oscillations the three modes settle 
into their quasi-stationary states, which change only sec- 
ularly as the spin and temperature of the star evolve. 
The modes heat the star viscously at constant Cl in seg- 
ment a —>■ b of the trajectory for ihcat-up ~ 1,100 yr. 
At point b, the star reaches a state of thermal balance. 
In leg b — > c the star continues its evolution in thermal 
equilibrium and slowly spins up due to accretion until 
the angular velocity evolution also reaches an equilib- 
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FIG. 3: (a)The amplitudes of the r-mode \C a \ and of the 
n = 13, m — —3 and n = 14, m = 1 inertial modes \Cp\ and 
C 7 | are shown as a function of time for a star that executes 
a cyclic evolution (same parameters as in Fig. [2}. The low- 
est parametric instability threshold is also displayed. (b)The 
ratio of the r-mode amplitude to the parametric instability 
threshold is plotted as a function of time. It can be seen that 
once the r-mode crosses the parametric instability threshold 
it remains close to it for the rest of the evolution. (c)The 
parametric instability threshold is displayed as a function of 
time. Its value changes as the angular velocity and tempera- 
ture evolve. 
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FIG. 4: The trajectory of a neutron star in the — T% phase 
space is shown for a model with T c = 5.0 x 10 9 K, / du = 0.03 
and Sns = 0.03 that reaches an equilibrium steady state. The 
star spins up until it crosses the r-mode stability curve and 
the r-mode becomes unstable. The r-mode then quickly grows 
to the first parametric instability threshold and excites the 
daughter modes. In leg a — > b of the trajectory the star is 
viscously heated by the mode triplet until the system reaches 
thermal equilibrium. Segment b — * c shows the star contin- 
uing to heat and spin up in thermal equilibrium until the 
accretion torque is balanced by the gravitational radiation 
emission. The r-mode stability curve represents the points 
in phase space where the viscous driving rate is equal to the 
gravitational driving rate. The H=C curve is the locus of 
points where the viscous dissipation due to the mode triplet 
balances the neutrino cooling. 
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FIG. 5: The (Q, Tg) initial values (region delimited by the 
solid line) that lead to equilibrium steady states and their 
corresponding final steady state values (region enclosed by 
the dashed line) are shown. Since both the initial and final 
values of T% are low, these evolutions are roughly independent 
ofT c . 



rium. The timescale to reach an equilibrium steady state 
is Steady ~ 3.5 x 10 6 yr for this set of parameters. 

Fig. [5] displays the possible initial values for the angu- 
lar velocity fl and temperature Tg of the star that lead to 
a balancing between the accreted angular momentum and 
the angular momentum emitted in gravitational waves. 
The fraction of the star that is above the threshold for 
direct URCA reactions and the slippage factor are varied 
within the corresponding "STEADY STATE" region of 
Fig. []Ib) . The final equilibrium values are also displayed 
and cluster in a narrower region than the initial values. 
Because viscosity is so small in this regime, the values 
of £1 also tend to be small. Thus, although an interest- 
ing physical regime, this case is most likely not relevant 
to recycling by accretion to create pulsars with spin fre- 
quencies as large as 716 Hz. Note that a steady state 
can be achieved when S ns = 0. This is the probable end 
state of the problem first calculated by Levin [i8| . The 
reason we do not find a cycle at low S ns is twofold: (1) 
the shear viscosity we are using is lower (shear viscosity 
in Ref. [l8| is amplified by a factor of 244), and (2) the 
nonlinear couplings keep all mode amplitudes small. 



C. Thermal Runaway Evolutions 

We now consider evolutions in which the three-mode 
system is not sufficient to halt the thermal runaway. We 
observe two such scenarios. In the first scenario, the 
star is unable to reach thermal equilibrium. The run- 
away occurs on a period much shorter than the accretion 
timescale and so the whole evolution is at approximately 
constant angular frequency. In the second scenario, the 
star reaches a state of thermal equilibrium but the spin 
evolution does not reach a steady state. The star contin- 
ues to spin up by accretion until it climbs to the peak of 
the H = C curve, thermal equilibrium fails and a run- 
away occurs. 



1. Fast Runaway 

A typical trajectory of a star that goes through a rapid 
thermal runaway is displayed in Fig. [5] This star has 
S ns = 0.25 and /du = 0.058. Initially, the growth of the 
r-mode is halted by the two daughter modes once the 
lowest parametric instability threshold is crossed, and 
the three modes settle in the (f2,T)-dependent quasi- 
stationary states of Eqs. ©. They viscously heat up the 
star until hyperon bulk viscosity becomes important for 
the daughter modes. As the amplitudes of the daughter 
modes decrease the coupling is no longer strong enough 
to drain enough energy to stop the growth of the r-mode. 
The daughter modes are completely damped and the r- 
mode increases exponentially. The system goes back to 
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the one-mode evolution described by Ref. 

Fig.[BJa) and (b) compare both the temperature evolu- 
tion and the trajectory in the fi— Tg, plane of the star for a 
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FIG. 6: This plot compares the full evolution resulting from 
solving Eqs. (|4|l. (fT3jl .rtT8 |l with the reduced Q. — T evolution 
that assumes the amplitudes go through a series of steady 
states Eqs. <[ig ]l -pP )l for a model with T c = 5.0 x 10 9 K, 
/du = 0.058 and S n s = 0.25. (a) The temperature is dis- 
played as a function of time for the two different methods, 
(b) The angular velocity Q = Q/£l c is shown as a function 
of temperature. The evolution occurs at constant spin fre- 
quency. It can be seen that the steady-state amplitude ap- 
proximation is extremely good. The 'X' shows the point at 
which the r-mode crosses its second lowest parametric insta- 
bility threshold, where additional dissipation would become 
operative. 



simulation solving the full set of equations to a simulation 
that assumes quasi-stationary solutions for the three am- 
plitudes and evolves only the angular velocity and tem- 
perature of the star. It can be seen that the steady state 
approximation is very good until the thermal runaway 
occurs. Afterward, the temperature evolution of the re- 
duced equations is offset slightly from the quasi-steady 
result and intersects the r-mode instability curve sooner. 
This evolution is similar to that described by Nayyar and 
Owen [24(. However, the r-mode crosses its second low- 



0.13 




FIG. 7: The trajectory of a neutron star in the — Tg 
phase space is shown for a model with T c — 5.0 x 10 9 K, 
/du = 4.0 x 10 and Sns = 0.02 that goes through a slow 
thermogravitational runaway. Portion a — > b of the trajectory 
shows the mode triplet heating up the neutron star through 
boundary layer and shear viscosity until the system reaches 
thermal equilibrium. Segment b — > c represents the accre- 
tional spin-up of the star in thermal equilibrium. The dotted- 
dashed line is the locus of points where the viscous dissipation 
of the mode triplet is equal to the neutrino cooling, and is la- 
beled as the H = C curve. The star reaches the maximum of 
this curve and fails to reach an equilibrium between the ac- 
cretion torque and gravitational emission. It then continues 
heating at constant angular velocity and crosses its second 
lowest parametric instability threshold, at which point more 
modes would need to be included to make the evolution accu- 
rate. Eventually the star reaches the r-mode stability curve 
again. 



est parametric instability much earlier in the evolution 
(see the 'X' in the figure), and at that point more modes 
need to be included to model the instability accurately. 
Thus, we cannot be sure that a runaway must occur in 
this case. We shall return to this issue in a subsequent 
paper. 



2. Slow Runaway 

In this section we examine evolutions in which the neu- 
tron star has both a very small slippage factor, S ns < 
0.03, and only a small percentage of the star is above the 
threshold for direct URCA reactions, fav < 5 x 10~ 5 . 
A trajectory for this kind of evolution is displayed in 
Fig. [71 After the star crosses the r-mode stability curve, 
the r-mode increases beyond the first parametric insta- 
bility threshold, and its growth is temporarily stopped 
by energy transfer to the daughter modes. As in the 
previous scenarios, the star is viscously heated by the 
mode triplet at constant Q in part a — > b of the trajec- 
tory on a timescale of about 5, 000 yr. At point 6, it 
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FIG. 8: The spin-down timescale is shown as slippage fac- 
tor S ns and fraction of the star subject to direct URCA fdv 
for cyclic evolutions are varied for a fixed hyperon critical 
temperature of T c — 5.0 x 10 9 K. This timescale dominates 
the heat-up timescale and hence represents the time the star 
spends above the r-mode instability curve. It increases as the 
viscosity is lowered and the star gets closer to the steady state 
region. 



reaches thermal equilibrium. In leg b — > c of the tra- 
jectory, the star continues its evolution by ascending the 
H = C curve and spinning up because of accretion for 
about 2 x 10 6 yr without finding an equilibrium state 
for the angular momentum evolution. Once it reaches 
the peak of the H = C curve, the cooling is no longer 
sufficient to stop the temperature from increasing expo- 
nentially and a thermal runaway occurs. The cross mark 
'X' on the trajectory shows the point at which the r-mode 
amplitude crosses its second lowest parametric instabil- 
ity threshold. At this stage more inertial modes need to 
be included to model the rest of this evolution correctly. 
As for the cases that evolve to steady states, these long- 
timescale runaways tend to occur at low spin rates. 
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FIG. 9: (a)The stability regions are plotted at fixed hyperon 
superfluidity temperature T c — 6.5 x 10 9 K, while varying /du 
and Sns- The steady state region remains roughly the same 
as in Fig. [ljb), the slow run-away region disappears, and the 
cycle region increases dramatically while shrinking the fast- 
runaway region, (b) The spin-down timescale is shown for the 
cyclic evolutions in part (a). 



V. PROBABILITY OF DETECTION 

Fig. [5] shows how the time the star spends above the 
r-mode stability curve changes when S ns and /du are var- 
ied. For large enough values of S na the boundary layer 
viscosity dominates. In this region of phase space the 
spin-down timescale can be approximated by 



t 



spin — down 



dt ~ 

—^dn 

b dfl 



(32) 



Jt£ r (4k) 2 o; / 3(D 7 1 AQ 

~<Q> 6 

in 



G 

250 yr 



AfkHz 



\Ca\ 2 
1 



<j/kHz> 7 Mi. 4 R| V |c 



where M XA = M/(1.4M ), R 6 = i?/(10 6 cm), 
VkHz — v/lkHz, I = 0.261 [17], the r-mode am- 
plitude at its parametric instability threshold |c^ h | « 

1.5 xlO- 5 , and 0, = Vfi|cJ 



\5u\/(Ak^/l^l^) 



14 



This approximation agrees with spin-up timescales ob- 
tained from our simulations to ~ 25%. 

The maximum v is approximately the same as the ini- 
tial frequency, and can be determined by equating the 
driving and damping rate of the r-mode, since it is on 
the r-mode stability curve 



800Hz 



Sri 



Mi A R 6 



4/11 



1 



2/11 



(33) 



Thus, the spin-down timescale is very sensitive to 
the slippage factor i spin -down oc 5 , ns 24/11 (A^ kHz /t'kHz)- 
The dependences on /du and accretion rate M are 
much weaker; a rough approximation, obtained by 
matching direct URCA cooling and nuclear heat- 



-1/9 
.4 I 

The 

wave amplitude measured at distance d [46 



ing, is T 8in cx M^f^R-^Ml 
sTfTM-^R-^Ml 



and 



34/99 



"max 

ravitational 
H is 



(34) 



Taking v « j/ max gives 

ho.S^M-^Rlfl^M-^K (35) 

The maximum distance at which sources could be de- 
tected by Advanced LIGO interferometers, assuming 

Kta = lO" 27 , H is 

d max « 3Mpc Mi. 4 Rl^ H z fej) (36) 

« 1.5 Mpc S^M^ 11 ^ 1711 



x T, 



-6/11 / Fa 



„th 



Eqs. (|33[) and (|36p imply that gravitational radiation 
from the r-mode instability may only be detectable for 
sources in the Local Group of galaxies. Eq. (|33f implies 
that for accretion to be able to spin up neutron stars to 
v > 700 Hz, we must require (S ns / 1 M lA R & ./%~) i / 11 > 1. 
Assuming this to be true, d max ^ 1-1.5 Mpc. However, 
^spin-down ~ 1000 yr at most, making detection unlikely 
for any given source. Moreover, unless S na can differ sub- 
stantially from one neutron star to another, only those 
with v given by Eq. (|33|) can be r-mode unstable. Slower 
rotators, including almost all LMXBs, are still in their 
stable spin-up phases. 

Still more seriously, Fig. \V[b) shows that spin cycles 
are only possible for S ns < 0.50, assuming T c w 5.0 x 10 9 
K; Eq. Q then implies v < 450 Hz. This would restrict 
detectable gravitational radiation to galactic sources, al- 
though the duration of the unstable phase could be 
longer. 



Within the context of our three mode calculation, 
S^s > 0.50, which is needed for explaining the fastest pul- 
sars, would imply fast runaway. There are two possible 
resolutions to this problem. One is that including addi- 
tional modes prevents the runaway; we shall investigate 
this in subsequent papers. The second is that T c is larger, 
or that neutron stars do not contain hyperons (e.g., be- 
cause they are insufficiently dense). Fig. [U(a) shows the 
same phase plane as Fig.QJb) but with T c = 6.5 x 10 9 K, 
and Fig. [5{b) shows the results for i S pin-down analogous 
to Fig. [5J Larger T c permits spin cycles for higher values 
of 5ns (and hence u), but the time spent in the unstable 
regime is shorter. 



VI. CONCLUSIONS 

In this paper, we model the nonlinear saturation of un- 
stable r-modes of accreting neutron stars using the triplet 
of modes formed from the n = 3, m — 2 r-mode and the 
the first two near resonant modes that become unstable 
{n = 13, m = —3 and n = 14, m = 1) by coupling to 
the r-mode. This is the first treatment of the spin and 
thermal evolution including the nonlinear saturation of 
the r-mode instability to provide a physical cutoff by en- 
ergy transfer to other modes in the system. The model 
includes neutrino cooling and shear, boundary layer and 
hyperon bulk viscosity. We allow for some uncertainties 
in neutron star physics that is not yet understood by 
varying the superfluid transition temperature, the slip- 
page factor that regulates the boundary layer viscosity, 
and the fraction of the star that is above the density 
threshold for direct URCA reactions. In all our evolu- 
tions we find that the mode amplitudes quickly settle 
into a series of quasi-stationary states that can be calcu- 
lated algebraically, and depend weakly on angular veloc- 
ity and temperature. The evolution continues along these 
sequences of quasi-steady states as long as the r-mode is 
in the unstable regime. The spin and temperature of 
the neutron star can follow several possible trajectories 
depending on interior physics. The first part of the evo- 
lution is the same for all types of trajectories: the star 
viscously heats up at constant angular velocity. 

If thermal equilibrium is reached, we find several pos- 
sible scenarios. The star may follow a cyclic evolution, 
and spin down and cool in thermal equilibrium until the 
r-mode enters the stable regime. It subsequently cools 
at constant f2 until it reaches the initial temperature. 
At this point the star starts spinning up by accretion 
until the r-mode becomes unstable again and the cycle 
is repeated. The time the star spends in the unstable 
regime is found to vary between a few hundred years 
(large S ns - 1) and 10 6 yr (small S ns - 0.05). Our 
cycles are different from those previously found by Ref. 
[18j in that our amplitudes remain small, ~ 10~ 5 , which 
slows the viscous heating and causes the star to spend 
more time in the regime where the r-mode instability is 
active. Furthermore, we find that the star stops heating 
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when it reaches thermal equilibrium and not when the r- 
mode reaches a maximum value. The cycles we find are 
narrow with the spin frequency of the star changing less 
than 10% even in the case of high spin rates ~ 750 Hz. 
Other possible trajectories are an evolution toward a full 
steady state in which the accretion torque balances the 
gravitational radiation emission, and a very slow thermo- 
gravitational runaway on a timescale of ~ 10 6 yr. These 
scenarios occur for very low viscosity (S ns < 0.04). Al- 
though theoretically interesting, they do not allow for 
very fast rotators of ~ 700 Hz. 

Alternatively, if the star does not reach thermal equi- 
librium, we find that it continues heating up at constant 
spin frequency until it enters a regime in which the r- 
mode is no longer unstable. This evolution is similar 
to that predicted by Nayyar and Owen [24|. However, 
the r-mode grows above its second parametric instability 
threshold fairly early in its evolution and at this point 
more inertial modes should be excited and the three- 
mode model becomes insufficient. Modeling this scenario 
accurately is subject of future work. 

We have focused on cases with T c >5x 10 9 K. These 
are cases for which the nonlinear effects are substantial. 
In this regime, hyperon bulk viscosity is not important 
except for thermal runaways where we expect other mode 
couplings, ignored here, to play important roles. Fast ro- 
tation re quir es large dissipation, as has long been recog- 
nized [l8ll30l| and these models can only achieve v > 700 
Hz if boundary layer viscosity is very large. Alterna- 
tively, at lower T c < 3 x 10 9 K, large rotation rates can 
be achieved at r-mode amplitudes below the first para- 
metric instability threshold [l[. Nayyar and Owen found 
that increasing the mass of the star for the same equation 
of state makes the hyperon bulk viscosity become impor- 
tant at lower temperatures [24|. Conceivably, there are 
accreting neutron stars with relatively low masses that 
have lower central densities and small hyperon popula- 
tions. These could evolve as detailed here and only spin 
up to modest frequencies. Hyperons could be more im- 
portant in more massive neutron stars leading to larger 
spin rates and very small steady state r-mode amplitude 
as found by Wagoner [l|. 

Our models imply small r-mode amplitudes of ~ 10~ 5 
and therefore gravitational radiation detectable by ad- 
vanced LIGO interferometers only in the local group 
of galaxies up to a distance of a few Mpc. The r- 
mode instability puts a fairly stringent limit on the 
spin frequencies of accreting neutron stars of u max w 
800Hz[S ns /(Mi A R e )]^ n Tg 2/n . In order to allow for 
fast rotators of > 700 Hz in our models a large bound- 
ary layer viscosity with (Sns/^i^-^VTsin) 4 / 11 ~ 1 is 
required. Slippage factors of order ~ 1 lead to time peri- 
ods on which the r-mode is unstable with a timescale of at 
most 1000 yr, which is about 10~ 3 times shorter than the 
accretion timescale. This would mean that only about 1 
in 1000 LMXBs in the galaxy are possible LIGO sources. 
However, lower slippage factors lead to a longer duration 
of the gravitational wave emission, but also lower fre- 



quencies. We also note that in this model we have con- 
sidered only very fast accretors with M ~ 10 _8 MQyr _1 
and most LMXBs in our galaxy accrete at slower rates. 
Investigations with more accurate nuclear heating models 
are a subject for future work. 

Our analysis could be made more realistic in several 
ways, such as by including the effects of magnetic fields, 
compressibility, multi- fluid composition [48[, superfluid- 
ity, superconductivity, etc. These features would render 
the model more realistic, but its generic features ought 
to persist, since the upshot would still be a dense set of 
mode frequencies exhibiting three mode resonances and 
parametric instabilities with low threshold amplitudes. 
Although the behavior of the star would differ quanti- 
tatively in a model different from ours in detail, we ex- 
pect the qualitative behaviors we have found to be ro- 
bust, as they are well described by quasi-stationary mode 
evolutions whose slow variations are determined by com- 
petitions between dissipation and neutrino cooling, and 
accretion spin-up and gravitational radiation spin-down. 
In our model, it seems that three mode evolution involv- 
ing interactions of the r-mode with two daughters at the 
lowest parametric instability threshold is often sufficient 
to quench the instability. Our treatment is inadequate to 
follow what happens when the system runs away; for this, 
coupling to additional modes is essential. For this regime, 
a generalization of the work of Brink et al. (2(| [27], [28[ 
that includes accretion spin-up, viscous heating and neu- 
trino cooling would be needed. Such a calculation is 
formidable even in a "simple" model involving coupled 
inertial modes of an incompressible star. 
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APPENDIX A 

This appendix will sketch the derivation of Eqs. (Q]) 
from the Lagrangian density. We follow closely Appendix 
A in Schcnk et al., which contains the derivation of the 
equations of motion for constant Q. 

The Lagrangian density as given by Eq. (Al) in Schcnk 
et al. [H is 

£=^e+^-B-£-^-C-£ + a 0Xt (i)-£, (A _i) 
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where the operators B • £ = 21~i x £ and 
p(C-Oi = -V J (r 1 pV^ J )+V lP V^ J +pV i ^(A-2) 

with rot = — (1/2) (S~2 x x) 2 . We are interested in a 
situation where the uniform angular velocity of the star 
changes slowly on the timescale of the rotation period 
itself. In order to remove the time dependence we define 
the new displacement and time variables 



dr = fldt. 



(A-3) 



In terms of these new variables the Lagrangian density 
can be written as 



C 



1 



a cxt (t) Z 



where the primes denote derivatives with respect to t, 
B = Vl~ 1 B and C = £1~ 2 C. The momentum canonically 
conjugate to £ is 



7T 



dC 



i'+Cixi 



(A-5) 



The associated Hamiltonian density is 

2 



H 



1. 



7T 



-B £ 



2^ 141 



2^ s 



a ext 

(A- 



Hamilton's equations of motions can be written as 



C' = T-C + F(r), 



(A-7) 



where 



7T 

the operator T is T = Tq + T± with 



To = 



±B 2 C 



1 



and 



and 



Ti = 



F(r) 







(VP)" 
\/fi 







n 3 / 2 



We assume solutions of the form £(t, x) = e iwt £(x). Spe- 
cializing to the case of no forcing term a ex t = leads to 
the eigenvalue equation 



(T - «D)C(x) = 0. 



(A-8) 



Since the operator To is not Hermitian it will have dis- 
tinct right and left eigenvectors. Similar to Schenk et 
al. [2!| we label the right eigenvectors of T as (a, and 
the associated eigenfrequencies as Co a = wa/O, and the 
eigenvalue equation above becomes 



(T - mDa)Ca(x) = 0. 
The left eigenvectors \A satisfy 

(Tj - iQ A )x A = 0, 



(A-9) 



(A-10) 



where 



T T 

J- n 



±B ±B 2 C 



1 



For simplicity, in this appendix we specialize to the case 
of no Jordan chains when the set of right eigenvectors 
forms a complete basis. The orthonormality relation be- 
tween right and left eigenvectors is 



XaXbJ = J d 6 xp(x)x A • Cb = 5 
We can expand £(r, x) in this basis as 
C(r,x)=^CU(r)a(x) 



AB ■ 



(A-ll) 



(A-12) 



where the coefficients C a are given by the inverse of this 
mode expansion 



Ca(t) = (xa,C(t-,x) 



(A-13) 



Using Eqs. (|B-2IA-9IA-11I) in Eq. (fFf)) leads to the equa- 
tions of motion for the mode amplitudes 



C A -iu A C A = 9(t)J2 C *b 



B \ XA, 







(A-14) 



+ (XA,F(r)), 



where g[t) = (VU)"/Vn. Following Sec. IV of Schenk 
et al. |29| we replace the externally applied acceleration 
by the nonlinear acceleration given by Eq. (4.2) of Ref. 
[29j. The inner product can be written in terms of the 
displacement variable £. The left eigenvectors are 



XA 



OA 
TA 



where ta can be chosen to be proportional to £a because 
they satisfy the same matrix equation. 



ta = -ilA/h 



(A-15) 



which corresponds to Eq. (A-45) in Schenk et al. \2i 
with the proportionality constant Ida = 0,~ 1 bA - 
MR 2 /Co A (also given by Eq. (2.36) of Ref. 0). 
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The equations of motion for the mode amplitudes be- 
come 



C A — iluaCa 



iMR 2 
b A 



J2 C b d 3 xi* A -i B (A-W) 



BC 



where the nonlinear coupling kabc — kabc I (MR 2 fl 2 ) 
and kabc is explicitly give by Eq. (4.20) of Ref. [29]. The 
g(r) integral mixes only modes with ttia = vn B because 
of the e lm ^ dependence of the displacement eigenvectors 
£. (f** d<j)e^ mA - mB ^ = if m A + m B .) So, this term 
will be zero for our mode triplet. Also, in the case of a 
single mode triplet there is only one coupling and Eqs. 
(|A-16[) take the form of Eqs. dJ). 



APPENDIX B 

In this appendix we study the behavior of the mode 
amplitudes and temperature near equilibrium assuming 
constant angular velocity. We are performing a first order 
expansion of Eqs. ([5]) and (fT8j) . Similar to Ref. [43], each 
of the five variables is expanded about its equilibrium 
(Xj) e as follows 

Xj(f) = {\C a \, \C P \, |C 7 U,T 8 } = (Xj) e [l + 0(f)] 

(B-l) 

where the perturbation \Q\ << 1 and j = a, /3, 7, T. The 
expansion leads to a first order differential equation for 
each Q 



df 



df 



df 



df 



df 



(a — C@ — C7 



\7a 
Cl\SCj 

(5 

V7/3 

(77 )e 

n 



tan ( 



(B-2) 



1 



Ca - C/3 + C7 + 

(^) e CT ' 
Ca + Cp - C7 + 

dT 8 

Ca 



tan i 



Ct 

la + 7/3 



tan i 



■7 7 



Co 



Co 



+ c- 



:tan0 e \ Q\Su\ 

"7a +7/3 - 7 7 \ , (7a 



C/3- 



"7a -7/3+77 



7/3 -7 7 ) 



n\6u>\ 
MR 2 n 2 c w,3% 



2k 2 WaWpUl-y I SQ I C(T e )T§ e 

2 ( W a -^Ca + W/3C/3 + W7C7 
V 7a 



n\5Q\ 
1 



tan i 



la dT 8 
fdL v \ I 



1 97 Q _ 1 97/3 
h — — 

7/3 <??8 

Ct, 



1 d7 7 

' 77 5T 8 



Ct 



n c fJ|£w|C(T e ) 

where the equilibrium amplitudes \Cj\ e have been written 
in terms of the corresponding driving and damping rates 
using Eqs. ([6]). Eq. (|B-2j) can be written in matrix form 
as 



ctr 

Let Cj oc exp(Af). The determinant \\A l - 
leads to the eigenvalue equation 

A 5 + a 4 A 4 + a 3 A 3 + a 2 A 2 + a x A + a = 0. 

The coefficients aj with j = 0, 4 are 

7/3 + 77 - 7a 



0.4 



"3 



«2 



«1 



2 tan< 
2 



n|fu| 
7§ + 7 2 + 7a 



(B-3) 

xs^\\ = 

(B-4) 
(B-5) 



tan 4> 2 


(f2|to|) 2 


la 7/377 


/ 12 


(fi|fc|) 3 


(tan^ + 


47a 7/37 7 


( 1 , 



tan ( 



-1, 



tan i 



"o 



(ft|<fo|) 3 \tsm<t> t 

2MR 2 tt 2 c ( 7 a7/37 7 ) 2 1 
K 2 u) a u)/3UjjC(T e ) (0|(5tD|) 4 tan</> e 

/ ^la\ tip f djp \ 
7a \dTsJ e 7/3 V^S/e 77 
47a 7/3 77 1 1 / 

(f2|to|) 3 tanc/> e n|<Jw|C7(T e ) V dT 8 
The eigenvalues can be approximated as 



1 



tan i 



Al,2 
A3, 4 

A 5 



«4 

~y ~ 

e ± iw, 
ai 



a 1 



a 4 
T 



(B-6) 



where e = (a 2 — et3a4)/a4 and w = yax/a^. The system 
is unstable when a 2 — 0304 > or ao < 0. The first 
two eigenvalues will have a negative real part as long as 
7/3 + 7y > 7a ■ If the heating compensates the cooling of 
the star ao ~ and becomes negative if the star can not 
reach thermal equilibrium. The other critical stability 
condition a 2 — 0304 = can be written as 



7a 



[i+r /J +r 7 -(r^+i^)-(r /J -r 7 ) J (r /J +r 7 )] = o, 

(B-7) 

where Tp = 7/3/7Q and T 7 = 7 7 /7 a . Note that we have 
ignored the smaller terms of order 0([7 Q /(f2|(5(I>|)] 5 ) 



This 
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condition can be rewritten by defining variables D\ = 
+ T 7 and D 2 — T p ~ T 1 

2 + 2L>i - D{ - D\ - 2D 2 2 D 1 = 0. (B-8) 

If D 2 = then the equation has one solution D\ = 1 + V3 
for D\ > 2, which corresponds to T = Tp = T 7 = 1.37 



and matches the result of Wersinger et al. [HI . For the 
viscosity we consider (see Sec. IIIDj) a 2 — a^a^ < 0. 
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